Disclaimer 1: The present material is an adaptation of the original R script prepared by Dr. Matteo Fontana for the a.y. 2020/2021 Nonparametric statistics course, and later Prof. Andrea Cappozzo (a.y 2022-2023) While I acknowledge Matteo and Andrea for the (great) work done I hereby assume responsibility for any error that may be present in this document.

Disclaimer 2: I will start from the assumption that you are all intermediate R users, if this is not the case please let me know, and we shall find a solution together.

Loading necessary libraries

library(MASS)
library(rgl)
library(DepthProc)
library(hexbin)
library(aplpack)
library(robustbase)
library(MDBED)  # plot exponential bivariate distribution

Computing depths with R: some general ideas

Let us start by simulating \(500\) bivariate datapoints whose distribution is marginally exponential

set.seed(2781991) # reproducibility
n=500
df_bivariate_exp = cbind(rexp(n), rexp(n))
head(df_bivariate_exp)
##           [,1]       [,2]
## [1,] 0.3004277 0.35894610
## [2,] 2.4430839 0.22392999
## [3,] 1.8444128 0.15064919
## [4,] 0.2383688 1.12290176
## [5,] 1.9613656 0.07966151
## [6,] 5.2480885 0.01762207

and visualize their scatterplot

plot(df_bivariate_exp[,1],df_bivariate_exp[,2], xlab="exp 1", ylab="exp 2")

We can further employ a hexagonal binning plot to visualize the data density

bin=hexbin(df_bivariate_exp[,1],df_bivariate_exp[,2], xbins=10, xlab="exp 1", ylab="exp 2")
plot(bin, main="Hexagonal Binning")

Depth measures (multivariate)

Now on depths: there are many possible packages hosted on CRAN to work with depths in a multivariate setting. We will hereafter use DepthProc, which is a good general purpose package and in addition it has good plotting capabilities. Even though not all depth measures introduced in class are directly available (simplicial, Oja, convex hull peeling depth), it includes the two depth measures we will be using throughout this section, namely Tukey and Mahalanobis depths. Let us look at the help file for the depth function.

Calculating depth for a given dataset is immediately accomplished by typing

set.seed(20)
tukey_depth=depth(u=df_bivariate_exp,method='Tukey')

It may be useful (spoiler alert, it will become apparent why when we study onparametric forecasting) to calculate the depth of a point relative to a sample. You can do it by:

depth(u = c(0.5, 0.5), X = df_bivariate_exp, method = 'Tukey') 
## Depth method:  Tukey 
## [1] 0.254

Compute the median (deepest point) with depthMedian function

depthMedian(df_bivariate_exp,depth_params = list(method='Tukey'))
## [1] 0.8364350 0.9297927

Or, having already computed the Tukey depth for the entire sample:

df_bivariate_exp[which.max(tukey_depth),]
## [1] 0.8364350 0.9297927

If you have the luck of dealing with a bivariate dataset, you can easily visualize the depth surface in a very convenient way. DepthProc offers you two possible methods:

DepthProc::depthContour(df_bivariate_exp,depth_params = list(method='Tukey'))

or

depthPersp(df_bivariate_exp,depth_params = list(method='Tukey'))

For additional special effects:

depthPersp(df_bivariate_exp,depth_params = list(method='Tukey'),plot_method = 'rgl')
Computational caveat

Since the Tukey depth is obtained by counting points in the space of all hyperplanes(surfaces) in \(\mathbb{R}^d\), when \(d>2\), the computational burden increases. Libraries usually recur to heuristic methods which lead to poor estimates (we will see this in the functional data depth measures lab)

a = depthPersp(df_bivariate_exp,depth_params = list(method='Tukey'))

The very same analysis can be carried out by considering a different depth measures:

Simplicial (or Liu) depth:

mrfDepth::sdepth(df_bivariate_exp)$depthZ
##   [1] 0.112827245 0.043805925 0.052095903 0.099119154 0.025286042 0.006000000
##   [7] 0.038261101 0.146898085 0.008824782 0.189377599 0.141395176 0.015819446
##  [13] 0.015610691 0.230319965 0.006000000 0.144052394 0.235806746 0.154630176
##  [19] 0.178300939 0.144672139 0.041322742 0.052061617 0.006516310 0.006000000
##  [25] 0.064843373 0.035115774 0.151183041 0.032581259 0.168724533 0.075996040
##  [31] 0.076487240 0.014556776 0.017281648 0.215428978 0.077749620 0.087816307
##  [37] 0.228942318 0.149266485 0.231516141 0.228164812 0.021211387 0.091886037
##  [43] 0.174609122 0.031402033 0.009396093 0.137307482 0.047747785 0.013995847
##  [49] 0.064876500 0.136080933 0.141716493 0.202258493 0.038173793 0.162751189
##  [55] 0.167804525 0.156815511 0.022217688 0.011608470 0.022862399 0.038947340
##  [61] 0.125680228 0.071467224 0.048179395 0.150906536 0.024388391 0.009300529
##  [67] 0.028471207 0.006373760 0.065631649 0.010840766 0.010340343 0.086454065
##  [73] 0.087399039 0.018473574 0.084852645 0.029433952 0.014314750 0.051906077
##  [79] 0.040759398 0.142870174 0.087026970 0.112785668 0.200215032 0.025810271
##  [85] 0.013132723 0.220927832 0.053925152 0.062973948 0.092953039 0.019528310
##  [91] 0.175128136 0.020190453 0.083810416 0.247232054 0.027455151 0.093649757
##  [97] 0.131388898 0.037729483 0.113098534 0.058061472 0.032074945 0.007998358
## [103] 0.034047179 0.017662892 0.006189294 0.026123234 0.029044692 0.127519231
## [109] 0.088771422 0.043270927 0.121903132 0.227215733 0.025504455 0.015007026
## [115] 0.217532318 0.185444914 0.010363619 0.036335273 0.068984234 0.029277398
## [121] 0.020245551 0.194301760 0.043730256 0.068852162 0.036250042 0.153333366
## [127] 0.174292633 0.021759905 0.177407490 0.141731077 0.144373663 0.157634208
## [133] 0.105141367 0.112218992 0.079795446 0.046371007 0.036009561 0.154733032
## [139] 0.128767994 0.030082816 0.075373735 0.046179540 0.009522611 0.065790859
## [145] 0.215487553 0.131210904 0.108150276 0.014290364 0.125446314 0.142174614
## [151] 0.091941908 0.237619287 0.234031871 0.149322162 0.210289398 0.134839317
## [157] 0.009364995 0.042450829 0.112996692 0.097948910 0.076924548 0.159482531
## [163] 0.009424729 0.020996547 0.032876355 0.014878432 0.094244730 0.043231137
## [169] 0.018420842 0.009918681 0.175285945 0.006648043 0.006000000 0.037928242
## [175] 0.192059010 0.096792621 0.028786295 0.006823671 0.082096192 0.196624671
## [181] 0.018690586 0.203099645 0.062382983 0.130567448 0.203177584 0.089904097
## [187] 0.063008620 0.013599826 0.047570466 0.025388609 0.132268151 0.223749620
## [193] 0.010616124 0.181931574 0.146335997 0.188846705 0.026875921 0.112860565
## [199] 0.008519931 0.166297656 0.075174445 0.010163315 0.050812130 0.238632784
## [205] 0.008371683 0.032559094 0.089776227 0.007795060 0.188754569 0.006000000
## [211] 0.219730401 0.012509211 0.193274887 0.078380810 0.132731777 0.091073859
## [217] 0.164242895 0.009730835 0.226219089 0.067238622 0.090459859 0.015671681
## [223] 0.019475143 0.007183620 0.036713910 0.148235217 0.052316682 0.022229664
## [229] 0.149962383 0.045797040 0.166040611 0.006071710 0.008082623 0.100052491
## [235] 0.087995557 0.234032064 0.047979960 0.006095420 0.011830746 0.066110776
## [241] 0.131828138 0.013506483 0.038720139 0.156290460 0.007411981 0.014224111
## [247] 0.056601685 0.228979549 0.065634884 0.211382089 0.182183065 0.080841104
## [253] 0.029272811 0.032994712 0.043482435 0.035512181 0.024525919 0.018416302
## [259] 0.212896057 0.035479827 0.043720356 0.246222276 0.056309921 0.157564720
## [265] 0.246196200 0.234541951 0.150499408 0.063904387 0.225730449 0.023263249
## [271] 0.010770505 0.010604293 0.121982568 0.032744380 0.068786151 0.112709564
## [277] 0.126749209 0.131819350 0.038298235 0.029107227 0.172471642 0.144846367
## [283] 0.011919115 0.035989425 0.059580172 0.018198518 0.016025545 0.244591786
## [289] 0.123340512 0.227583504 0.014720477 0.090419586 0.124363909 0.033742231
## [295] 0.169045851 0.009196658 0.173329985 0.009098969 0.172314074 0.038714972
## [301] 0.082954149 0.243680325 0.222993747 0.006000000 0.019619818 0.101167298
## [307] 0.206802328 0.011461574 0.011817563 0.022908274 0.130585412 0.045820895
## [313] 0.006000000 0.042339040 0.072620228 0.171963590 0.012823913 0.217507159
## [319] 0.131091436 0.119595722 0.063843011 0.028823430 0.015158268 0.211571094
## [325] 0.006634764 0.245146631 0.222121061 0.114795326 0.023690176 0.168820050
## [331] 0.175471425 0.190391723 0.030776831 0.088089142 0.246577058 0.119702731
## [337] 0.017343844 0.202930970 0.176901369 0.132935027 0.010365985 0.009306855
## [343] 0.010799334 0.212724051 0.162729797 0.026289446 0.023675930 0.038878383
## [349] 0.013983485 0.050880991 0.170724727 0.031582104 0.065803607 0.078078036
## [355] 0.112293116 0.067284304 0.171121424 0.006000000 0.006000000 0.063771447
## [361] 0.228675906 0.193270686 0.198568993 0.132339329 0.108548326 0.110324746
## [367] 0.057557959 0.113764010 0.124622691 0.009608760 0.016955212 0.016568511
## [373] 0.008966173 0.171838569 0.009098100 0.225907091 0.017107516 0.020021006
## [379] 0.075964169 0.021267933 0.083915349 0.093347080 0.217454330 0.204529348
## [385] 0.102622353 0.008511433 0.184174421 0.006213970 0.204944202 0.201729193
## [391] 0.096672381 0.006000000 0.084770312 0.097965425 0.042973224 0.044008209
## [397] 0.073059082 0.026855977 0.220022696 0.248786537 0.071409856 0.197489775
## [403] 0.109334718 0.070016756 0.058649830 0.047266002 0.032400753 0.081144892
## [409] 0.021029674 0.204261197 0.078311418 0.007849337 0.006000000 0.096135162
## [415] 0.011256682 0.175142816 0.101600116 0.083376874 0.065253640 0.224461212
## [421] 0.006000000 0.011608567 0.038358017 0.032649299 0.077289132 0.057689403
## [427] 0.013298549 0.140689234 0.015727069 0.006402685 0.174051911 0.241760147
## [433] 0.006000000 0.205702924 0.148022358 0.014803197 0.015973489 0.048148152
## [439] 0.025040153 0.088875389 0.073833933 0.026484342 0.015265905 0.204336432
## [445] 0.037277302 0.204826279 0.103561822 0.041502861 0.241818094 0.225796267
## [451] 0.097488664 0.006047903 0.020918222 0.036242075 0.094477533 0.047743197
## [457] 0.149372915 0.058438419 0.074095082 0.065148466 0.008410218 0.031138035
## [463] 0.194908854 0.203296038 0.032386315 0.129884444 0.173873772 0.060589661
## [469] 0.029493010 0.111200135 0.139690658 0.184330299 0.057656808 0.007309800
## [475] 0.228035203 0.168880315 0.055882898 0.055419755 0.020802666 0.144387281
## [481] 0.028750320 0.080912669 0.243760485 0.107519666 0.085550281 0.122074704
## [487] 0.138308907 0.071798006 0.118166212 0.032619214 0.037605235 0.106827148
## [493] 0.087675206 0.122298766 0.025643190 0.068792766 0.092982688 0.007552551
## [499] 0.102314364 0.009563851

and now Mahalanobis depth

maha_depth <- depth(df_bivariate_exp,method='Mahalanobis')

This can be easily hard coded, recall the definition you have seen in class:

\[{M}_{h}D({F;X} ^ {n}) = \frac{ 1 }{ 1 + {{(x - \mu_F)} ^ {T}}{{\Sigma} ^ {-1}}(x - \mu_F) }\]

sample_mean <- colMeans(df_bivariate_exp)
sample_S <- cov(df_bivariate_exp)

maha_depth_manual <- 1/(1+mahalanobis(x = df_bivariate_exp,center = sample_mean,cov = sample_S))

And check that the obtained result is equal to the one obtained via the depth function ( beware of comparing floats)

all(abs(maha_depth-maha_depth_manual)<1e-15) # food for thought: sqrt(2) ^ 2 == 2?
## [1] TRUE
all.equal(as.vector(maha_depth)+1, as.vector(maha_depth_manual))
## [1] "Mean relative difference: 0.6672736"
depthMedian(df_bivariate_exp,depth_params = list(method='Mahalanobis'))
## [1] 0.9211997 1.0228290
df_bivariate_exp[which.max(maha_depth_manual),]
## [1] 0.9211997 1.0228290

And again the graphical outputs:

depthContour(df_bivariate_exp,depth_params = list(method='Mahalanobis'))

depthPersp(df_bivariate_exp,depth_params = list(method='Mahalanobis'))

depthPersp(df_bivariate_exp,depth_params = list(method='Mahalanobis'),plot_method = 'rgl')

Please note that anything that comes out of depthContour or depthPersp is NOT a density: we are not doing density estimation here, we are performing data exploration by means of a nonparametric procedure. Let us plot the probability density function (NOT the depth function) of the bivariate exponential distribution we drew a sample from:

PDF_3dPlot(rho = 0, Betax = 1, Betay = 1, title="PDF (NOT depth function) of the bivariate exponential distribution", zlabel="")

depthPersp(cbind(rexp(1e4), rexp(1e4)),depth_params = list(method='Tukey'), plot_title="Tukey depth for the exp distribution")

We notice that the tukey depth concentrates around the median, where as the exponential bivariate function increases as x,y go to 0.

Exercise

Try out the routines seen so far on something even more exotic:

set.seed(1992)
df_bivariate_cauchy = cbind(rcauchy(n,location=0,scale=.001), rcauchy(n,location = 0,scale=.001))
head(df_bivariate_cauchy)
##               [,1]         [,2]
## [1,] -0.0024184800 0.0008798105
## [2,] -0.0006282804 0.0009932694
## [3,] -0.0006975300 0.0010347342
## [4,] -0.0038080204 0.0141419885
## [5,]  0.0061098712 0.0014929490
## [6,] -0.0012696290 0.0008865863

The Cauchy distribution is a distribution that has neither the first moment (the mean) nor the second moment (the variance). This means that the CLT does not apply.

  • Can we still perform non-parametric data exploration with depth measures?
  • Would the Mahalanobis depth be a sensible measure to be used in this context? Why not?

Multivariate outlier detection via depth measures

In the previous Section we have computed depth measures for a given dataset and we have seen how to effectively plot them. This resulted in a convenient way to nonparametrically explore a (bivariate) dataset. Nevertheless, one of the main aims for a statistician to employ depth measures is to perform multivariate outlier detection. To appreciate this, let us simulate \(100\) data points, of which \(95\%\) comes from a multivariate normal with mean vector mu_good and covariance matrix sigma_common and the other 5% from another multivariate normal, which we assume is our outlier generator process.

mu_good = c(0,0) 
mu_outliers = c(7,7)

sigma_common = matrix(c(1,.7,.7,1), ncol = 2)

frac = .05
n=100
# sample points
n_good=ceiling(n*(1-frac))
n_outliers=n-n_good
df_contaminated_normals = data.frame(rbind(
  mvrnorm(n_good, mu_good, sigma_common),
  mvrnorm(n_outliers, mu_outliers, sigma_common)
))

Let us visualize the true nature of our dataset

class <- c(rep(1,n_good),rep(2,n_outliers))
plot(df_contaminated_normals,xlab="Norm 1", ylab="Norm 2",col=class)

We can clearly see those red points up there… how can we flag them in an automated fashion? The depth contour plot surely helps

depthContour(
  df_contaminated_normals,
  depth_params = list(method = 'Tukey'),
  points = TRUE,
  colors = colorRampPalette(c('white', 'navy')),
  levels = 10,
  pdmedian = F,
  graph_params = list(cex=.01, pch=1),
  pmean = F
)

But, as we have seen in class, a very handy graphical tool for spotting multivariate outliers is the bagplot. The bagplot function from the aplpack package can be used to display bagplots for bivariate data

bagplot(df_contaminated_normals, factor = 3, show.whiskers = T)

Remember:

aplpack::bagplot(df_contaminated_normals,show.whiskers = F,main="Bagplot")

aplpack::bagplot(df_contaminated_normals,show.loophull = F,main="Sunburst plot")

In addition, if we save the output of the bagplot to an object, we can automatically extract the outliers

bagplot_cont_normals <- bagplot(df_contaminated_normals)

outlying_obs <- bagplot_cont_normals$pxy.outlier

Once the outlying units have been identified, one can discard them from the original data and keep working on the clean subset only. There are several ways to do this.

A more “sql-oriented” approach (this is just a code alternative):

df_clean_1 <-
  dplyr::anti_join(
    x = df_contaminated_normals,
    y = outlying_obs,
    by = c("X1" = "x", "X2" = "y"),
    copy = TRUE
  )

A more “object-oriented programming” approach:

ind_outlying_obs <- which(apply(df_contaminated_normals,1,function(x) all(x %in% outlying_obs)))
df_clean_2 <- df_contaminated_normals[-ind_outlying_obs,]
all.equal(df_clean_1,df_clean_2)
## [1] TRUE

Outlier detection in Star Cluster CYG OB1 data

Data for the Hertzsprung-Russell Diagram of the Star Cluster CYG OB1, which contains \(47\) stars in the direction of Cygnus, from C.Doom.

data(starsCYG, package = "robustbase")
names(starsCYG)
## [1] "log.Te"    "log.light"
plot(starsCYG, main="Star Cluster CYG OB1")

We can see two groups of points: the majority which tends to follow a steep band, the so called Main Sequence, and four stars in the upper-left corner. In astronomy the \(43\) stars are said to lie on the Main sequence and the four remaining stars are the red giants, namely points with indexes 11, 20, 30 and 34. In details, the red giants are very bright, but they emit light with a very low color-temperature (and thus their surface temperature is still fairly low). We can easily isolate them thanks to the procedure seen so far. Let us look at the perspective plot of the depth surface

depthContour(as.matrix(starsCYG), depth_params = list(method='Tukey'), points=TRUE)

As you can see, the sample mean is biased due to the presence of the 4 red giants, whereas the Tukey Median is not. Let us compute it:

depthMedian(starsCYG) 
##    log.Te log.light 
##      4.38      5.02

As before, we can use the bagplot to visualize and flag the outlying data points.

bagplot_starsCYG <- with(starsCYG,aplpack::bagplot(log.Te,log.light))

red_giants <- bagplot_starsCYG$pxy.outlier
ind_outlying_obs <- which(apply(starsCYG,1,function(x) all(x %in% red_giants)))
ind_outlying_obs
## [1] 11 20 30 34

Moving beyond \(R ^2\)

So far, we have only dealt with bivariate data. Even though the computational complexity escalates quickly when it comes to compute depth measures in \(R ^d\), with \(d>2\), we can still appreciate their usefulness in moderate dimension. Let us generate a trivariate dataset with contamination.

mu_good = rep(0,3)
mu_outliers = c(12,12,3)

sigma_common = diag(3)*2

frac = .1
n=300
# sample points
n_good=ceiling(n*(1-frac))
n_outliers=n-n_good
df_3 = data.frame(rbind(
  mvrnorm(n_good, mu_good, sigma_common),
  mvrnorm(n_outliers, mu_outliers, sigma_common)
))
class <- c(rep(1,n_good),rep(2,n_outliers))
pairs(df_3, col=class)

To visualize the outliers in this context we retort to a bagplot matrix:

bagplot_matrix <- aplpack::bagplot.pairs(df_3)

You can notice that outlier detection becomes more difficult when the dimension increases, as

These phenomena are respectively denoted as masking and swamping, we will cover it in details during the robust statistics module of this course!

Let us conclude this lab session by looking at some DD-plots. For two probability distributions \(F\) and \(G\) , both in \(R ^d\), and \(D(\cdot)\) an affine-invariant depth, we can define depth vs. depth plot being very useful generalization of the one dimensional quantile-quantile plot:

\[ D D(F, G)=\left\{\left(D_{F}(x), D_{G}(x)\right) \text { for all } x \in \mathbb{R}^{d}\right\} \]

Its sample counterpart calculated for two samples \(\mathbf{X} =\left\{X_1, \cdots, X_n \right\}\) from \(F\), and \(\mathbf{Y} =\left\{Y_1, \cdots, Y_m \right\}\) from \(G\) is defined as

\[ D D\left(F_{n}, G_{m}\right)=\left\{\left(D_{F_{n}}(z), D_{G_{m}}(z)\right), z \in\{\mathbf{X} \cup \mathbf{Y}\}\right\} \]

The ddPlot function automatically computes and plots it:

df_good <- df_3[1:n_good,]
df_out <- df_3[(n_good+1):n,]
ddPlot(x = df_good,y = df_out,depth_params = list(method='Tukey'))
## DDPlot

## 
## Depth Metohod:
##   Tukey

They would line up if they came from the same distribution

It is easy to manually build a DD-plot

depth_good <- depth(u = df_3,X = df_good,method = "Tukey")
depth_out <- depth(u = df_3,X = df_out,method = "Tukey")
plot(depth_good,depth_out, col="blue", xlab="X depth", ylab="Y depth", main= "Depth vs. depth plot")
grid(10, 10, col="grey50", lty=1)
abline(0,1, col="grey50")

Clearly the distributions are not the same… What if we had some extra samples coming from \(F\) (e.g., the trivariate normal centered at \(0\))?

n_extra <- 100
df_extra <- data.frame(mvrnorm(n_extra, mu_good, sigma_common))
ddPlot(x = df_extra, df_good,depth_params=list(method='Tukey'))
## DDPlot

## 
## Depth Metohod:
##   Tukey

Indeed this time we can conclude the same distribution generated the data.

Measuring robustness.

Breakdown point

fracs <- seq(0.1, .5, length.out=25)
set.seed(41703192)
n_good <- 100
df.clean <- data.frame(mvrnorm(n_good, mu_good, sigma_common))
multivar.mean.clean <- colMeans(df.clean)
multivar.median.clean <- depthMedian(df.clean,depth_params = list(method='Tukey'))

mean.diffs <- numeric(length(fracs))
median.diffs <- numeric(length(fracs))

for (i in 1:length(fracs)){
  n=100
  # sample points
  frac <- fracs[i]
  n_good=ceiling(n*(1-frac))
  n_outliers=n-n_good
  
  df_contaminated_normals = data.frame(rbind(
  as.matrix(df.clean[1:n_good, ]),
  mvrnorm(n_outliers, mu_outliers, sigma_common)
  ))
  current.mean <- colMeans(df_contaminated_normals)
  current.median <- depthMedian(df_contaminated_normals, depth_params = list(method='Tukey'))

  mean.diffs[i] <- max(abs(current.mean - multivar.mean.clean) )
  median.diffs[i] <- max(abs(current.median - multivar.median.clean))
}

plot(x=fracs, y=mean.diffs, main="Empirical Breakdown point")
lines(x=fracs, y=median.diffs, col="red")

Empirical Influence function

sim.data <- rnorm(20)
stat1 <- function(data){base::mean(data)}
stat2 <- function(data){median(data)}

L = 25
grid.for.x <- seq(-1.5*min(sim.data), max(sim.data)*1.5, length.out=L)

stat1.0 <- base::mean(sim.data)
stat2.0 <- stat2(sim.data)

EIF1 <- numeric(L)
EIF2 <- numeric(L)
for (i in 1:length(grid.for.x)){
  sim.data.i <- sim.data
  sim.data.i[1] <- grid.for.x[i]
  EIF1[i] <- abs(base::mean(sim.data.i) - stat1.0)
  EIF2[i] <- abs(stat2(sim.data.i) - stat2.0)
}

plot(x=grid.for.x, y=EIF1, ylim=c(-.05, .12))
lines(x=grid.for.x, y=EIF2, col="red")

An alternative is the empirical influence value: the (estimated) the derivative of a sample estimator with respect to each statistical unit in the sample of its input. One must be careful on the definition of the statistic. It is implemented in the boot package and the statistic (i.e. a functional of the sample), has to be an (R) function of the data matrix and a vector of weights, i.e the weight associated to each statistical unit (row in the matrix), see more information here

For simplicity, we try it with the \(t\) statistic on univariate data: \(t(\mathbb{X}) = \frac{\hat{\mu}(\mathbb{X})-\mu_0}{\hat{\sigma}(\mathbb{X})}\) where \(\mathbb{X}\) denotes the data sample, with \(hat{\mu}(\mathbb{X})\) denoting the sample mean, \(\mu_0\) a hypothesised mean and \(\hat{\sigma}(\mathbb{X})\) the sample standard deviation

library(boot)
## 
## Caricamento pacchetto: 'boot'
## Il seguente oggetto è mascherato da 'package:robustbase':
## 
##     salinity
mean.statistic <- function(data, w){
  return ()
}